Locally critical point in an anisotropic Kondo lattice 
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We report the first numerical identification of a locally quantum critical point, at which the criti- 
cality of the local Kondo physics is embedded in that associated with a magnetic ordering. We are 
able to numerically access the quantum critical behavior by focusing on a Kondo-lattice model with 
Ising anisotropy. We also establish that the critical exponent for the q— dependent dynamical spin 
susceptibility is fractional and compares well with the experimental value for heavy fermions. 



How to properly describe heavy fermion metals near 
quantum critical points (QCPs) is a subject of intensive 
current research. It is well established experimentally [1] 
that these systems are prototypes of non-Fermi liquid 
metals [2,3]. In a number of cases, striking deviations 
from the commonly applied T — spin-density- wave pic- 
ture (usually referred to as the Hertz-Millis picture) [4] 
have been seen. In particular, experiments [5-9] have 
shown that the spin dynamics in the quantum critical 
regime can display fractional exponents essentially ev- 
erywhere in the Brillouin zone, as well as uj/T scaling. 
These features are completely unexpected in the Hertz- 
Millis picture - which corresponds to a Gaussian fixed 
point - and they directly imply the existence of quantum 
critical metals that have to be described by an interacting 
fixed point. A number of theoretical approaches are be- 
ing undertaken to search for such non-Gaussian quantum 
critical metals [10-13]. 

Here we are concerned with a new class of QCP [10,11] , 
which has properties that bear a close similarity to those 
seen experimentally. The key difference between the tra- 
ditional Hertz-Millis QCP and this locally critical point 
(LCP) is that, in the latter, the local Kondo physics itself 
becomes critical at the antiferromagnetic ordering tran- 
sition. Such a LCP was shown to arise in an extended 
dynamical mean field theory (EDMFT) of a Kondo lat- 
tice model. The latter was mapped onto a self-consistent 
impurity model - the Bose-Fermi Kondo model - which 
in turn was analyzed using a renormalization-group (RG) 
approach, based on an e expansion. 

One of the key issues is whether the destruction of the 
Kondo effect is accompanied by a fractional exponent in 
the frequency/temperature dependences of the dynami- 
cal spin susceptibility. The fractional nature of the expo- 
nent has been seen experimentally, and is characteristic 
of a non-Gaussian magnetic quantum critical metal. The 
calculation of this exponent requires a detailed numerical 
study of the quantum critical behavior. 

In this paper we demonstrate that, for an anisotropic 
version of the Kondo lattice model, the EDMFT equa- 
tions can be solved using a Quantum Monte Carlo 
(QMC) method [15], including at its QCP. We also ana- 
lyze them using a saddle-point approximation [16,17]. 



We consider the following Kondo lattice model: 
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Here, ty specify the bandstructure (cu) and density of 
states [po(e)] of the conduction c— electrons, Jk describes 
the on-site Kondo coupling between a spin-i local mo- 
ment Si and a conduction-electron spin density s Cj j and, 
finally, Iij denote the exchange interactions between the 
z— components of two local moments. 

We study this model using the EDMFT ap- 
proach [18-20], which maps the model onto a self- 
consistent anisotropic Bose-Fermi Kondo model: 
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where the parameters E p , w p and g are determined by a 
set of self-consistency equations. The latter are dictated 
by the translations! invariance: 
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where G\ oc (u>) is the local conduction-electron Green 
function and X1oc(t) =< T t S z (t)S z (0) > is the local 
spin susceptibility. The conduction-electron and spin 
self-energies are and M(ui), respectively, with 
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where the Weiss field Xo 1 characterizes the boson bath, 
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It has been recognized [10,11] that the solution depends 
crucially on the dimensionality of the spin fluctuations, 
which enters through the RKKY-density-of-states: 
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(6) There is also a constraint equation that reads 



In this paper, we will consider the following specific form, 
characteristic of two-dimensional fluctuations: 



Pl (e) = (1/27)6(7- |e|), 
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where is the Heaviside function. Enforcement of 
the self-consistency condition on the conduction-electron 
Green function is not essential for the critical properties 
discussed here, provided that the density of states at the 
chemical potential (/i) is finite; the corresponding bath 
density of states ^ p 5(oj — E p ) = N is also finite. 

If, instead of the general form Eq. (6) we were to choose 
a semi-circular RKKY-density-of-states (representative 
of the 3D case), our EDMFT equations would become 
essentially the same as those of Rcfs. [15-17,21]. 

To cast the impurity model in the form of a func- 
tional integral, we adopt a bosonized representation of 



the fcrmionic bath, cl = F\- 



and a canonical 

transformation using U = exp[— i(v / 2/2) ( f> s (S (T crX crcr )], 
where [22] X a p = \a >< f}\ are the Hubbard operators 
and $ s = (<I>| - $|)/\/2, leading to 
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Here, 7io(0s, <t>c) describes the bosonized conduction elec- 
tron bath, S* = (X u f/f t + H.c.)/2, a' = {*fc) x=0 h 
is the spin density of the conduction electron bath, T = 
^ and T z = 3$ [tan- (^) - f ] . (T and T z 
can vary independently when we allow the longitudinal 
and spin-flip parts of the Kondo interaction to be differ- 
ent.) Eq. (8) describes an Ising spin in a transverse field 
r, with a retarded self-interaction that is long-ranged in 
time [15]. The associated partition function is 
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Here, g c ~ (l/AT )[ln N 2 j ] 1 and JC C comes from inte- 



grating out the electron bath (k c ~ r^ATg) 
JC c (iu n ) = K c \w n \. 
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To first gain qualitative insights, we carry out a saddle- 
point analysis of Eq. (9). This analysis is formally exact 
when the number of components for the field n is gener- 
alized from 1 to AT and a large-A" limit is subsequently 
taken [16,17]. At the saddle point level, iX = A , and 
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The QCP is reached when the static susceptibility at the 
ordering wavevector Q, specified by 7q = —7, becomes 
divergent. Given that [18-20] 



X (q, W ) = [A7H+7 q ]- 1 , 
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it implies M(u> = 0) — > 7. This, in turn, establishes that 
Xioc(w = 0) is also divergent [through Eqs. (3,7)] and that 
Ao = Xo = I [from Eqs. (11,4)]. The self-consistency 
equation (3) then simplifies and becomes 
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= 7 + 27exp[-27xioc(w)]. (14) 
We immediately see, from Eqs. (11,12,14), that 
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are self-consistent to the leading order. 

Since the w-dcpendent part of M(u>) is subleading com- 
pared to both l/xioc(w) and the w-dependent part of 
Xq (u>), we can determine the spin self-energy entirely in 
terms of the leading order solution for xioc 

(uj) in Eq. (16) 

and the second equality in Eq. (14) : 



M(u>) = 7 + 27(-^/A) c 
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with the exponent a = I/A. Now, to satisfy the first 
equality of Eq. (14), we need to add sub-leading terms 
to the Weiss field and replace Eq. (15) by 



X^M =7-2Aln- 1 
Eq. (11) then leads to: 
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At the saddle point level, however, the critical ampli- 
tudes for the local susceptibility are pinned to the initial 
parameters of the Weiss field [23]. The self-consistent 
amplitudes of xioc is not determined from the leading 
terms alone and d = c; a fractional a cannot emerge. 
We now show that a fractional a does arise [23] in the 
physical case corresponding to Eq. (8). 

We have directly studied the physically relevant case 
[Eq. (8)] numerically using the QMC algorithm of 
Refs. [15,24]. Starting from a trial XoH 1 ")' we com- 
pute xioc(t). A new X(7 1 ( r ) is then obtained from the 



2 



FIG. 1. Static local susceptibility vs. T for I/T% = (+), 
0.75 (*), 1 (x), 1.0625 (□), 1.125 (O), 1.1875 (O) and 1.25 
(A). The solid lines are the fits described in the text. Inset: 
results for I = 1.1875 T%, with 8(T = 0) = 5 x 10~ 3 . 

self-consistency Eq. (3), and the process is repeated un- 
til convergence is achieved. The results reported below 
are obtained for n c = tt [cf. Eq. (10)]. In this special 
case, the model can be solved exactly at I = 0, pro- 
viding a check on the QMC algorithm [15]. We choose 
NqT = 0.19 such that the bare Kondo scale (the inverse 
of the static local susceptibility at / = 0), NqT% = 0.17, 
is well above the lowest temperature attainable in our 
simulations, T min = 10~ 2 T%. The number of Trotter 
time slices used varies between 64 (for the higher tem- 
peratures) and 512 (for the lowest one). We perform be- 
tween 10 4 and 10 6 MC steps per time slice and between 
five and twenty self-consistency iterations. The numeri- 
cal error is controlled by the last step. We estimate that 
the calculated susceptibilities are accurate to within 3%. 

Fig. 1 shows the temperature dependence of Xk>c(w — 
0) for several values of /. The results exhibit very little 
/-dependence above T w 0.1 T%. At lower temperatures, 
a saturation of xioc at an /-dependent value is seen in 
the five lower curves. The two upper curves show no 
saturation; instead, the results are consistent with a log- 
arithmic T-dependence as shown in the inset of Fig. 1, 
which contains data at many more points of temperature. 
To interpret these numerical results, and inspired by the 
arguments developed earlier, we use the following fitting 
function for Xu>c( u ) ] 

XLM = tJ© (A - M) \ tanh (^) , (20) 

where a, A, and 5 are fitting parameters that may depend 
on / and T. The parameter 5 is a measure of the proxim- 
ity to the QCP as Eq. (20) implies that x _1 (Q) = 2/(5" 
for 5 <C 1. By Hilbcrt transforming Eq. (20) we obtained 
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FIG. 2. Matsubara-frequency dependence of xioc at 
T = 0.011 T% and I = 1.1875 T%. The solid line is the fit ob- 
tained from the Hilbert transform of Eq. (20). Inset: the same 
data on a logarithmic plot with the point at ui n — excluded. 
The dashed line is the predicted critical behavior [10,11] with 
a = 0.72 and A = 1.54 Tg-. 

an analytic expression for xioc(^„) that we used to an- 
alyze the numerical results. 

In the parameter region T < 0.1 T%, 1 < I/T% < 1.25 
all our data can be fitted with a single value of a = 
0.72 ± 0.01 and A = 1.54 T%. The error bar represents 
the amplitude of variation of a when it is allowed to 
freely adjust for each value of / and T in the critical re- 
gion. The fits are of excellent quality as shown in Fig. 2 
where we display results obtained for / = 1.1875 T% at 
T = 0.011 T%. The fitted values of xiocK* = 0) are 
represented by the solid lines in Fig. 1. We found in ad- 
dition that, within this range of values of / and T, the 
fitted S(I,T) can be described by the phcnomcnological 
expression S cx [5 (I) + ^S'^(I) + 4(T/T%) 2 }/2 that can 
be derived using Eq. (20) in the normalization condition 
(12). The parameter Sq, that decreases linearly with in- 
creasing /, allows us to determine the location of the 
QCP from the criterion 5(I C ,T = 0) = 0: I c w 1.2 T%. 

We have also computed X _1 (Q) ; the inverse of the 
static peak susceptibility. Fig. 3 shows that x _1 (Q) oc 
^0.72 f or j an( j ji j n q uan t um critical regime. Our 

numerical value for a [26] is very close to that seen [5] 
in the Ising-like system CeCu6- x Au x . The region I > I c 
will be discussed elsewhere [25]. 

We stress that a simultaneous treatment of the Kondo 
and RKKY couplings is crucial for our conclusion that 
the LCP solution is indeed self-consistent. 
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FIG. 3. Scaling plot of x _1 (Q) near the QCP. 



Some questions have recently been raised concerning 
the matching of the logarithmic terms in the LCP solu- 
tion [11,14] in the isotropic model. Burdin et al. [14] 
carried out large-iV and numerical analyses for zero 
Kondo coupling. They found that the spin liquid (SL) 
phase [28,29] is unstable at low T and conjectured that 
the LCP solutions may not be self-consistent either. 
Since the SL phase corresponds to the stable fixed point 
of the effective impurity model [27] this analysis covers a 
different parameter regime and is complementary to ours. 
By directly accessing the unstable fixed point, our results 
establish that the logarithmic terms are self-consistent in 
the Ising case. Whether numerical and analytical (be- 
yond RG) studies in the isotropic case will yield a self- 
consistent LCP similar to what we have shown here for 
the anisotropic model is left for future work. 

In summary, we have numerically identified a LCP so- 
lution in an anisotropic Kondo lattice model. The ex- 
ponent for the q— dependent dynamical susceptibility is 
fractional and is close to the experimental value. 
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